A boundary element regularised Stokeslet method applied to 

cilia and flagella-driven flow* 



David J. Smitht 

School of Mathematics and School of Clinical and Experimental Medicine, 
University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom 
Centre for Human Reproductive Science, Birmingham Women's NHS Foundation Trust 
Metchley Park Road, Edgbaston, Birmingham B15 2TG, United Kingdom 



Abstract 



A boundary element implementation of the regularised Stokeslet method of Cortez is applied to 
cilia and flagella-driven flows in biology. Previously-published approaches implicitly combine the 
force discretisation and the numerical quadrature used to evaluate boundary integrals. By contrast, 
a boundary element method can be implemented by discretising the force using basis functions, and 
calculating integrals using accurate numerical or analytic integration. This substantially weakens 
the coupling of the mesh size for the force and the regularisation parameter, and greatly reduces the 
number of degrees of freedom required. When modelling a cilium or flagellum as a one-dimensional 
filament, the regularisation parameter can be considered a proxy for the body radius, as opposed 
to being a parameter used to minimise numerical errors. Modelling a patch of cilia, it is found 
that: (1) For a fixed number of cilia, reducing cilia spacing reduces transport. (2) For fixed patch 
dimension, increasing cilia number increases the transport, up to a plateau at 9 x 9 cilia. Modelling a 
choanoflagellate cell it is found that the presence of a lorica structure significantly affects transport 
and flow outside the lorica, but does not significantly alter the force experienced by the flagellum. 
Key words: regularised Stokeslet, cilia, flagella, boundary element, slender body theory 

1 Introduction 

Fluid dynamic phenomena on microscopic scales are fundamental to life, from the feeding and swim- 
ming of bacteria and single-celled eukaryotes to complex roles in the organ systems of mammals. 
Examples include sperm swimming, ovum and embryo transport, respiratory defence, auditory per- 
ception and embryonic symmetry-breaking. Cilia, which range from a few to tens of microns in length, 
are hair-like organelles project from cell surfaces either in dense mats or one-per-cell, and through spe- 
cialised side-to-side or rotational beating patterns and viscous interaction with the cell surface, propel 
fluid in a parallel direction. Flagella may be tens to hundreds of microns in length, typically occurring 
singly or two per cell, and generally produce propulsion tangential to the flagellum axis through the 
propagation of a bending wave, resulting in swimming as in sperm, or the generation of feeding cur- 
rents as in choanoflagellates. These phenomena have since the work of Taylor, and Gray and Hancock 
in the 1950s, motivated the study of methods for solving the zero Reynolds number linear 'Stokes 
flow' equations, 
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and 



= -Vp + fiV 2 u + f 
V • u = 0, 




where u, p, \x and / are velocity, pressure, kinematic viscosity and force per unit volume respectively. 
These equations model flow strongly dominated by viscous effects with negligible inertia, valid for 
cilia and flagella since the Reynolds number can be estimated as Re = pJJLjp < 10 -2 . Hancock 



(1953), working with Sir James Lighthill, introduced the singular 'Stokeslet' solution, corresponding 



to a concentrated point-force, or equivalently the purely viscous component of the flow driven by a 
translating sphere. The Stokeslet is the solution of the Stokes flow equations for unit force acting in 
the j-direction and concentrated at corresponding to taking fix) = 5(x — £)ej, where 5(x — £) 
denotes the Dirac delta distribution and ej is the appropriate basis vector. The i-component of the 
velocity field driven by this force is written as Sij(x,£,), and in an infinite fluid takes the form 



Sij{x, 40 



21 -L % 3 
3 



(2) 



where r 7 - 



■ & and r 2 = \x — £| 



+ rn + t%, with <5jj denoting Kronecker delta tensor. The flow 



due to a force F concentrated at the point £ corresponds to taking f(x) = 5(x — £)F in equation Q, 
the solution being given by Ui(x) = (l/8ir(J,)Sij(x,£)Fj, We assume the summation convention over 
repeated indices throughout. 

The Stokeslet forms the basis for slender body theory for Stokes flow (Burgers, 1938; Hancock, 1953 



Tuck 



1964 



Batchelor, 1970 



example Pozrikidis 



Johnson, 1980), its local simplification 'resistive force theory', introduced 



1992 



by Gray and Hancock (1955) and the single-layer boundary integral method for Stokes flow (see for 

These techniques have proved very useful in the modelling of cilia- 



2002) 



driven flows (Blake, 1972, Liron and Mochon, 1976; Gueron and Liron, 1992, Smith et al. , 2007), sperm 



motility (Dresdner et al. 1980; Smith et al. 2009b), bacterial swimming (Higdon 1979b 



Phan-Thien 



et al. 1987 Ramia et al. 1993) and atomic force microscopy ( Clarke et al.[ |2006[ ), allowing three- 
dimensional flow problems with moving boundaries to be solved with greatly reduced computational 
cost compared with direct numerical solutions of the differential equations ([!]). These techniques 
involve expressing the fluid velocity field with an integral equation of the form 



«(*)~//(0-S(*,s-)d% 



(3) 



where S is a collection of surfaces or lines, for example representing boundaries, cell surfaces, cilia or 
flagella, with the symbol J s denoting line or surface integrals as appropriate, and /(£) denoting force 
per unit area or length respectively. The force exerted by the body on the fluid is given by /(£)d<f?£, 
and the force exerted by the fluid on the body is given by— f(£)dS£. The force density is calculated 
by combining equation ^ with either prescribed motions of the bodies, or a coupled elastic/active 
force model, through the no-slip, no-penetration condition = £• 



2 The method of regularised Stokeslets 

Line distributions of Stokeslets, for example equation Q with the boundary S parametrised as £(s) for 
scaled arclength parameter < s < 1, representing objects such as cilia or flagella, have the principal 
disadvantage that the flow field is singular at any point x = To calculate the force per unit 

length, it is necessary to perform collocation at points on the 'surface' of the slender body, displaced a 
small distance from the centreline X(s g ) = £( s <?) + a(s q )n(s q ), where a(s q ) is equivalent to the slender 
body radius and n(s q ) is a unit normal vector. The velocity field is regular in the region designated 
the outside of the body, with the singular line distribution lying inside the notional surface of the 
cilium. 

Point distributions of Stokeslets ^2 q= iF q ■ S(x,£ q ), for example representing a swarm of cells, 
are also singular at any x = £ . Surface distributions of Stokeslets, as used in the single-layer 
boundary integral method for Stokes flow, do not result in singular velocities, but nevertheless require 



2 



careful numerical implementation in order to ensure that the surface integrals are calculated correctly 



(Pozrikidis, 1992 2002) 



In order to circumvent these difficulties and ensure a regular flow field throughout the flow domain, 



Cortez (2001) introduced the 'regularised Stokeslet'. This is defined as the exact solution to the Stokes 



flow equations with smoothed point-forces, 



and 





V • u 



-Vp + nV 2 u + fip e (x-£), 



0. 



(4) 



The symbol ip t denotes a cutoff- function or 'blob' with regularisation parameter e, satisfying j R3 ip e (x)dV x 



Cortez et al. (2005) showed that with the choice ip t (x— £) := 15e 4 / (87r/irJ), the regularised Stokeslet 

5ij(r 2 + 2e 2 ) + nr 



velocity tensor is given by 



(5) 



Here and in the rest of the paper, we use the compact notation r € = yr + e . We follow Cortez 
and co-authors in using the solution Sfj, and its counterpart near a no-slip boundary given in 



equation (12), as the basis for our computational study. 



Cortez et al. (2005) derived the equivalent Lorentz reciprocal relation, and hence a boundary 



integral equation for the RSM. This leads to the following equation for the fluid velocity at location 
x, where the surfaces and lines representing the cells and beating appendages are denoted S, 



u[x 



-L//(e.s-(*,«)d*. 



(6) 



As above, the unknown /(£) is a force per unit area or length depending on whether £ is on a surface 
or line comprising S. A significant advantage of the RSM is that the kernel remains regular for x = £, 
even when S consists of lines or points rather than surfaces, and so collocation can be performed 
at such points without the need to make a finite displacement. This method has proved effective 
in modelling a number of biological flow problems, including flow due to an individual cilium and 



the swarming of bacteria (Ainley et al. 2008), the bundling of bacterial flagella (Flores et al. 2005 



Cisneros et al. 2008 ) , hydrodynamic interaction of swimming cells ( Cisneros et al. , 2007 ) and the 



flagellar motility of human sperm (Gillies et al. 2009). 



Alternative implementation of the RSM as a boundary element 
method 



The mathematical problem considered in this paper is the calculation of the force density /(£) from 
known boundary velocity, given by the time derivative through applying the no-slip condition 
u (£) = £ to equation Q. The boundary, representing a cell surface, a flagellum, an array of cilia, 
a network of fibres or a collection of particles, is discretised as a set of points {£ q : q = 1, . . . ,N}. 
Applying collocation at these points gives the following equation 



(7) 



Cortez et al. 


(2005 


); 


Ainley et al. 



(2008); Gillies et al. (2009); Cisneros et al. (2007) is equivalent to approximating the force density by 



constant values, /(£ r ) 
giving 



/ r , and replacing the integral of the Stokeslet with a low order quadrature, 



1 N 

— Y 

r=l 



(8) 
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where {wi, . . . ,wn} are quadrature weights. Hence the matrix equation AX = B can be formed, 
where 



^3(<z-l)+i,3(r-l)+j = u; r-5'| j (^,^ r )/(87r/i), 
-^3(r-l)+j = (fr)j 

and B 3{q _ 1)+i = u»(£ 5 ). (9) 

This discretisation can be viewed as a constant-force implementation of a boundary integral 
method, utilising a low order quadrature in which the abscissae are identified with the collocation 
points. However, while the force density /(£) varies relatively slowly, the kernel S e (£ g ,£) will vary 
rapidly in the neighbourhood of £„. Hence if the surface is discretised as patches Si, . . . , SV, the term 
w r f r ■ S e (£ g ,£ r ) in equation Q may not be an accurate approximation to f s /(£) ■ S e (£ 9 ,£) dS^, 
particularly when r = q. 

Consequently, in order to obtain accurate solutions, it is necessary to use a relatively large number 
of nodes £„, many degrees of freedom, and hence a large dense matrix system in order to calculate 
accurate solutions, as discussed in detail in [7|p2]). This is an unnecessary computational effort, since 
even with a constant-force discretisation, far fewer degrees of freedom for / are required than quadra- 
ture nodes for the evaluation of the kernel integrals. It is also necessary to employ an empirical rule of 
the form e = C(5s) m relating the regularisation parameter to the mesh size 5s to obtain the expected 



Cortez et al. ( 


2005 


); Ainley et al. 


(2008) 



For this reason, we instead implement a constant-force boundary element method where the inte- 
gration and force discretisation are 'decoupled', which we then apply to a number of example problems. 
On each patch S r , the force is approximated by f r . Then we have 

1 N r 

u ^ = E fr ■ / S£ (^> o d % (io) 

This problem can again be written as a matrix equation AX = B, with A now taking the form 

A 3(q-l)+i,3(r-l)+j = <Sf/(£ 9 ,£) d <% (H) 

The e-5s coupling of the original method is now substantially weakened, for the following reason: the 
discretisation chosen for the force no longer has any bearing on the evaluation of the kernel integrals 



in equation (11). If numerical integration is used it is still necessary to choose a sufficiently refined 
method for the value of e chosen, for the 'near-singular' integrals. However, for the values of e typically 
used, this is far less computationally costly than increasing the number of degrees of freedom, as shown 
in [7| and moreoever may be accomplished in principle with an automatic integration routine or in 
some cases analytic integration. 

The choice of constant-force elements is purely for simplicity: higher-order methods may also be 
used, generalising the above so that f(s) is approximated as ^2^ = i 4> r (s)f r , where . . . , <?W(s)} is 



a set of basis functions. Cubic Lagrange interpolating polynomials were used by Smith et al. (2009a) in 
the context of the non-regularised viscoelastic Stokeslet distributions. The integrals J s S e (£ g ,£) d5^ 
may be evaluated using any appropriate means: Gauss-Legendre numerical quadrature, more advanced 
adaptive quadratures, or in the case of problems where the boundary is reduced to a line segment, the 
integrals may be performed analytically (see [7]). In [7] we compare the two approaches ^ and (10). 



4 The Ainley et al. image system for a no-slip plane boundary 

In order to model fluid propulsion by cilia protruding from a cell surface, |Blake| ( |1971[ ) derived the 
Stokeslet image system satisfying the no-slip, no-penetration condition u = on X3 = 0. The technique 
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(a) cilium (2,2), timestep 37/40 

Fluid velocity component, 



(b) cilium (3, 1), timestep 11/40 

Fluid velocity component, «i(£(s)) 



Collocation point velocity, f(s s ) ■ Collocation point velocity, £(s q ) 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

arclength s arclength s 



Figure 1: Post-hoc verification of the accuracy of the solution, comparing the fluid velocity on the 
cilium tt(£(s,t n )) for dimensionless arclength < s < 1 and the prescribed cilium velocity on the 
collocation points £{s q ,t n ) for q = 1,...,N, at two different timesteps t n . Results were computed 
with a square lattice of 3 x 3 cilia with spacing 0.08L, beating in synchrony, for one beat cycle, with 
N = 20 collocation points and Nt = 40 timesteps per cycle. Results are shown for (a) an 'average' 
case, cilium (2, 2), timestep 37/40, for which the absolute RMS error is 0.0132, and (b) the worst case 
in the beat cycle, on cilium (3,1), timestep 11/40, for which the error is 0.0357. All velocities are 
given in dimensionless units, scaled with respect to oL. 



is equally useful in modelling the interaction of a glass surface and a sperm cell, and removes the need 
for the cell surface boundary to be included in the boundary integral equation. We shall not repeat 
the solution here, however it consists of a Stokeslet in the fluid, an equal and opposite image Stokeslet, 
a Stokes-dipole and potential dipole. As an important step in generalising the regularised Stokeslet 



method to the modelling of biological flows, Ainley et al. (2008) derived the equivalent regularised 



image system, which we denote Bf- for flow near a no-slip boundary. We rewrite the solution in index 



notation for direct comparison with the result of Blake (1971): 



S lj (r 2 + 2e 2 ) + 



TiT 



i' J 



+2hA jk 
6he 2 



S ij (R 2 + 2e 2 ) + R i R j 



d fhRi 5 i3 (R 2 + 2e 2 ) + RiR 3 



dR k V #f 

(SaRj-SijRa). 



4mh5 ik 4> e {R) 



(12) 



The tensor A 7 fc takes value +1 for j = k = 1, 2, value —1 for j = k = 3 and zero otherwise, originally 



w ritten as (5j a 5 a k — Sj3$3k) by Blake. The last line of equation (12) is not precisely equivalent to that 
Ainley et al. (2008), the latter having a typographical sign error. The term (j)e(R) := 3e 2 /(47ri?^) is 



m 



a more slowly-decaying blob than ip e (R), discovered by Ainley et al. (2008) as the function generating 



the potential dipole appropriate to the no-slip image system. The additional term on the final line is 
the difference between rotlets arising from the two different blobs ipe(R) and (p t (R), and is 0(e 2 /R 4 ). 



5 



(a), timestep 10/40 
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_1 -1 Qxi 1 

(c), timestep 30/40 
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(b), timestep 20/40 
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(d), timestep 40/40 

0.2 0.4 0.6 0.8 1 1.2 
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0.2 0.4 0.6 0.8 1 1.2 




Figure 2: Velocity field results for a square patch of 121 cilia, with (a-d) showing timesteps 10/40, 
20/40, 30/40 and 40/40 over one complete beat cycle. The upper panels show a 'vertical' section 
of the velocity field at X2 = 0.5; the lower panels show a 'horizontal' section at X3 = 0.5. Velocity 
magnitude is shown in colour, direction is shown with arrows. 
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(a), fixed cilia number 3x3, variable cilia spacing 




1.5 



0.5 1 

cilia spacing / L 

(b), fixed patch dimension 1.0L x l.OL, variable cilia number 



T3 

s 
= 




cilia spacing / L 

Figure 3: Mean volume flow rate Q versus cilia spacing, (a) for a patch of 3 x 3 cilia, the patch size 
varying from 0.16L x 0.16L to 3.0L x 3.0L, (b) for a patch of fixed size l.OL x l.OL, the cilia numbers 
varying from 2 x 2 to 13 x 13. 



5 Application to biological flows 



The numerical schemes given by equations ^ and ( 11 ) are applied to test problems involving rod and 
sphere translation in an infinite fluid in [7| In this section we apply the method to two biologically- 
inspired problems: determining the flow generated by an array of cilia protruding from a flat surface, 
and by a choanoflagellate cell in an infinite fluid surrounded by a silica lorica structure. We nondi- 
mensionalise with scalings L, 1/cr and aL for length, time and velocity, where L is cilium or flagellum 
length and a is radian beat frequency. The appropriate scaling for force per unit length on a filament 
is jiaL, and for force per unit area on a surface is jia. 



5.1 Modelling a patch of cilia 

Cilia lining the airway epithelia and the surfaces of microorganisms typically appear in dense 'fields', 
while cilia within the female reproductive tract generally appear in patches, with ciliated cells being in- 
terspersed between secretory cells. The density of cilia in such systems can be very high, approximately 



6-8/ /im 2 (Sleigh et al. 



1988), which for a square lattice is equivalent to a spacing of 0.35-0.41 /mi, or 

6 /an. This contrasts with embryonic nodal primary 



approximately 0.06L-0.07L for cilium length L 
cilia, for which there is only one cilium per cell. 

Slender body theory does not in general produce exact solutions, except in certain special cases 



such as that of a straight ellipsoids in a uniform flow, or uniform shear (Chwang and Wu, 1975). 



When modelling single bodies with non-uniform velocity or curvature, or multiple bodies, the flow 
field u(x,t n ) generated by a numerical solution for f(s,t n ) typically exhibits errors when compar- 
ing u(£(s,t n )) with the prescribed boundary movement £(s,t n ) for values of s that were not used 
as collocation points. These errors can be reduced in magnitude through mesh refinement, utilising 



higher-order singularities (Smith et al. 2007), or higher-order basis functions (Smith et al. 2009a), 



however they cannot be eliminated entirely at the ends of the body without resorting to replacing the 
line distribution by a surface distribution, which incurs considerably higher computational expense. 



In an earlier computational study utilising singular Stokeslet distributions, Smith et al. (2007) found 
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that cilia could not be arranged in arrays with a spacing of less than 0.3L, equivalent to approximately 
1.8 /im, without significant numerical errors occurring in the computational solution — a principal diffi- 
culty being that collocation points on the cilia may be close together, particularly during the recovery 
stroke. This limited the applicability of the technique in determining the effect of cilia density at 
realistic concentrations. 

In this study, we test the RSM without the use of higher-order singularities, and without the use 
of surface distributions, in order to model a densely-packed patch of cilia. The array is modelled as a 
square lattice of M = M\ x Mi cilia, for 1 ^ I ^ Mi and 1 ^ m ^ M% at time t the velocity field is 
given by, 



u(x, t) 



-. Ml M 2 ,-\ 

Mi M 2 N f . r / N 

- E E E # m (*) • / t l,m ( s > *)) ^, 



(13) 



where the constant vector /J.' m (i) denotes the force on the rth segment of the cilium in the (l,m) 
position in the array, at time t. 



We use the mathematical formulation of the tracheal cilium beat cycle given by Fulford and Blake 



(1986), based on the micrograph observations of Sanderson and Sleigh (1981). The mathematical 



specification has the following form: 



£) = E a p( s ) cos (pt) + P P ( S ) sin(pt), 
p=i 



(14) 



where the functions a p (s) and f3 p (s) are cubic polynomials, the coefficients being given in Fulford 
and Blake (1986). Rather than impose a metachronal wave on the cilia, we specify that they beat in 



synchrony, since fully-developed metachronal waves are generally more typical of larger cilia arrays. 
Hence the position vector of the cilium in the (I, m)-position is simply 



^ m (s, t) = t) + (l- l)aei + (m - l)ac 2 , 



(15) 



where the parameter o is the cilia spacing. 

In the numerical implementation, the flow field is calculated at t n = nSt, where n = 1, . . . ,Nt 
and 5t = 2-k/Nt- We find that results for the volume flow rate converge to three significant figures 
using N T = 40. As discussed in (7) and ff§A$, we choose the regularisation parameter based on the 
cilium radius to length ratio, so that e = 0.1/6 = 0.0167. Results are calculated using the University 
of Birmingham's BlueBEAR Opteron cluster, however the most computationally expensive simulation 
of 169 cilia requires only 8 hours 20 minutes of CPU time, and so could be performed on a desktop 
PC. 

To monitor the accuracy of the numerical solution for closely-approaching cilia, we use a similar 
technique to Gillies et al. (2009) and evaluate the post-hoc collocation error at each timestep t — t n , 



denoted £k m ■ This is given by 



?I,m 




)) 



t m (s,t n )\2ds, 



(16) 



where the integral is calculated numerically using the points s p = (p — l/2)/120 for 1 ^ p ^ 120. For 
all results presented, we used = 20 elements per cilium, and verified that the error did not exceed 
£n ,m < 0.04 on any cilium (/, to) or any timestep n. Example results showing the post-hoc comparison 
are given in figure [T] The discretised matrix A has similar favourable properties to the matrix resulting 
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from the discretisation using a singular kernel, with the diagonal entries being relatively large. The 
condition numbers calculated for the multi-cilia results were approximately 1.3 x 10~ 3 . 

Figure [2] shows examples of the three-dimensional flow field computed with the model at four points 
in the beat cycle. The rapid decay of the flow magnitude in both horizontal and vertical directions is 
evident. The tendency of a patch of cilia to draw fluid 'in' from each side to the rear, and push fluid 
'out' to the side in front is evident, as may be anticipated from the far- field form of the image system 



(Blake, 1971) 



An important function of cilia is the transport of liquid, and it is of interest to determine how 
the variation in cilia density observed in different biological systems might affect transport. In order 
to quantify this, we calculated the mean volume flow rate Q, using the formula described in [7| To 
examine the different fluid dynamic mechanisms at work, we simulate different cilia densities for flow 
driven by a finite patch in two ways: (a) keeping the number of cilia constant and varying the spacing 
through the size of the patch, (b) keeping the size of the patch constant, and varying the spacing 
through the number of cilia. Results are shown in figure [3} 

For a 3 x 3 cilia patch, as the cilia are brought closer together they hydro dynamically interact, 
reducing the fluid transport (figure [3^,). Because of this interaction effect, for a patch of fixed di- 
mensions, increasing the cilia number from 3 x 3 to 13 x 13 results only in a ~ 50% increase in fluid 
transport, and moreover 9x9 cilia give approximately the same transport as a 13 x 13 array (figure^). 
The common point in figures [3^,, b corresponds to cilia spacing 0.5L, for which the nondimensional 
mean volume flow rate is approximately 1.45. These results were computed with a prescribed beat 
pattern and frequency — in the biological system, these will vary depending on the flow field produced 
by the other cilia, as discussed in ^6} 



5.2 The modelling of a choanofiagellate 

Marine choanoflagellates are the nearest unicellular relatives of the animal kingdom, and exhibit a 
very diverse array of structures and behaviours that have yet to be fully explained biologically. Their 
feeding and swimming behaviour involves the beating of a single flagellum, and has been the subject of 



previous fluid mechanical modelling utilising slender body theory ( Orme et al. , 2003 ). For recent review 



of the biology, see Leadbeater et al. (2009), and for details on flagellar movement and fluid transport, 
see Pettitt et al. (2002). Certain species exhibit a silica surrounding basket-like structure, known as 



a lorica, made up of a network of silica strips, termed costae. The lorica has been hypothesised to 
perform a number of roles including protecting the cell, reducing the translational velocity of free- 
swimming cells, and directing feeding currents to its mouth. Simulating these behaviours is a subject 
for considerable future study; in this paper we simply demonstrate that with the boundary element 
regularised Stokeslet method it is computationally inexpensive to simulate such flows, and we make 
an initial study of the effect of surrounding the cell with a simple lorica structure. 

Both the flagellum and the 16 costae are modelled as filaments, the cell body being modelled as 
a sphere with 6x4x4 elements, with N = 30 constant-force elements for the flagellum, and N = 8 
constant-force elements for each costa. The costae were modelled less precisely than the flagellum, 
since their role is simply to provide a hydrodynamic drag. For these simulations we choose e = 0.01, 
corresponding to a flagellum aspect ratio of 0.1 /xm/10/im. The same regularisation parameter was 
also used for the sphere surface, since this has been verified to give accurate results, as shown in [7j 
and for the costae. Simulations for a single timestep were performed on a desktop PC, requiring under 
ten minutes of CPU time. 



Our simplified model for the choanofiagellate flagellum, based on that in Higdon (1979b), consists 
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of a sphere of radius a = 0.15, located at (0, 0, ho), where ho = 0.4. The flagellum is described by 



£l(s, t) = Ep(s) cos(ks — t), 
&(s,t) = E F (s) sin(«s - t), 

Us,t) = h + s, (0<s<l), (17) 

where t is scaled with respect to inverse radian frequency, and s is scaled with respect to flagellum 
extent rather than length. The wavenumber n = 2tt, so that the wavelength is unity. The enve- 
lope parameter E-p(s) = 1 — exp(— (k^s) 2 ), where ke = k. The coordinate s is not an arclength 
parametrisation, and so it is necessary to make a change of variable to arclength parameter s in order 
to calculate the force per unit length shown in figure [6j 

The silica lorica is constructed using a simplified model, based on the more complex family of 



geometric models of Leadbeater et al. (2009), in which it was shown that there are two sets of costae 



that run vertically or spiral around the axis of the cell, creating a basket-like structure. We defined a 
family of Nq vertical strips, given by £ V;m (s) for m = l,...N c , 

^' m (s) = E c (s)co S (2Trm/N c ), 
&' m {s) = E c (s)sm(27Tm/N c ), 

tl' m (s) = 2s, (0 < s < 1). (18) 
The family of Nq spiralling strips, denoted £ S ' m (s) for m = 1, . . . Nc, was then generated by changing 



the argument of the trigonometric functions in equation (18) to 2ir(m/Nc — s). For this study we 
chose the envelope function Ec(s) = 3.6s(l — s), and a total of 2Nc = 16 costae. 

Figure [4] shows the instantaneous fluid velocity field, computed on horizontal and vertical surfaces, 
both with and without the lorica. While the flow close to the flagellum is relatively unchanged by the 
presence of the lorica, it is apparent from figure |4]jc,d) that the remainder of the flow field is strongly 
suppressed in magnitude. Figure [5] shows the vertical component of the velocity field 113 only, on a 
horizontal cross-section. The presence of the lorica reduces this component in magnitude by nearly 
50% within the lorica interior, and again almost completely abolishes any external flow. By contrast, 
the hydrodynamic interaction of the flagellum and lorica is relatively weak, with figure [6] showing that 
the force density components are almost unchanged by its presence. 

6 Discussion and future work 

A simple modification of the regularised Stokeslet method was presented, based on decoupling the force 
discretisation and boundary integration. This reduced the number of degrees of freedom required to 
compute accurate solutions, and removed the need for an empirical parameter relating regularisation 
parameter and the discretisation mesh for the force. The implementation was designed to be the most 
elementary possible, to emphasise its ease of use in biological fluid mechanics problems. For the Stokes' 
law test example, a considerable reduction in computational cost was obtained, with approximately 
1000-fold fewer kernel evaluations being required to obtain accuracy to within 1%, and a 144-fold 
reduction in the number of degrees of freedom N. For problems involving M bodies, the storage and 
linear solver costs are at least 0(N 2 M 2 ), hence a reduction in the number of degrees of freedom for 
many-body problems will be very useful for the modelling of complex biological flows. 

The regularisation parameter e has two distinct roles in models that employ surface distributions 
and line distributions. For surface distributions, e plays the role of a tunable parameter, with smaller 
values giving more accurate solutions at the expense of more abscissae being required if numerical 
quadrature is used. It is sufficient to choose e small enough that the regularisation error is negligible; 



Cortez et al. (2005) showed for the blob function used in this paper that this error is 0(e 2 ), and in 



this study we found that for the test case of a translating sphere, a value of e = 0.01 gave errors of 
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Figure 5: Choanoflagellate (a) without, (b) with, a lorica structure consisting of 16 silica costae, 
showing the flow field U3 generated. 
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Figure 6: Choanoflagellate flagellar force components per unit length versus arclength s, with and 
without the 16 costae lorica structure. 
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less than 1% with 162 degrees of freedom. For line distributions, e is instead equivalent to a physical 
parameter of the filament — the slenderness ratio, as shown in[7j A detailed asymptotic analysis of the 
slender body theory shall follow this paper. 

The technique is capable of providing accurate results for denser arrays of cilia than have previously 
been possible, without recourse to higher order solutions or basis functions. Arrays of up to 169 cilia, 
with spacing 0.08 times their length could be simulated for a 40 timestep cycle in approximately 8 
hours with a single CPU. It was found that for prescribed cilia movement, cilia density of less than 
0.125 times their length does not produce significant increases in transport, due to hydrodynamic 
interaction obliterating any gains resulting from there being more cilia — it is likely that the role of 
realistically-high cilia densities will only be fully understood through more complex fluid-structure 
interaction models which take into account more biophysical details of the system in question. 
For discussion of the role of cilia in reproduction, 



sec 



Fauci and Dillon (2006). An additional 
complication when considering airway cilia is the effect of the serous-mucous bilayer, and possible 
pressure gradients which may exist due to surface tension, osmosis or both. Pressure gradients will 



have important effects on fluid transport (Mitran 2007, Smith et al. 2008b), and cilia penetration 



of the mucous layer may allow continued transport when many cilia are inactive (Fulford and Blake 



1986). 



The technique was also used to simulate the effect of the complex silica lorica structure of a marine 
choanoflagellate on the flow field and force density, with minimal computational cost. The lorica does 
not significantly affect the force distribution on the beating flagellum, but does significantly reduce the 
far-field flow produced. Future models which consider force and torque-balance for a free-swimming 
cell, in addition to ambient currents, and other features of the cell geometry such as the collar, may 
provide insight into the role and function of loricae, and the diversity that exists in choanoflagellate 
species. 

In line distribution representations, the force density depends on the regularisation parameter — in 
particular it can be estimated as ~ log(l/e) as e — > 0. For free-swimming problems, for example the 
work of Gillies et al. (2009), the swimming speed depends through force and torque balance on the 

this 



force density, and hence on the value chosen for the regularisation parameter. As shown in |7||1[ 
parameter may be chosen based on the slenderness ratio. This observation is relevant only to line 
distribution models of flagella, and not surface distributions. 

The calculations using line distributions were less accurate than those for translation of a sphere, 
however the method performed with similar accuracy to our previously-published studies for cilia 
motion, and superior accuracy for dense arrays. Future work may involve the use of higher-order 
solutions, higher-order basis functions and non-uniform meshes to reduce the end error for slender 
bodies. Recently, Smith et al. (2009a) presented a slender body theory for modelling linear viscoelastic 
flow, which exploits computational methods developed using singular solutions of the Stokes flow 
equations. A limitation is that errors in the slender body velocity field are propagated through the 
solution, which can lead to numerical instability if the Deborah number is larger than approximately 
0.06 for the problem considered. The RSM, possibly combined with higher order Green's functions, 
may provide a more accurate solution, particularly for problems with multiple cilia, and hence may 
extend the range of Deborah number that may be modelled. 

Other researchers have regularised the slender body theory equations differently: the 'Lighthill 
theorem' (Lighthill, 1976), later exploited by Dresdner et al. (1980) and by Gueron and Liron (1992); 



Gueron and Levit-Gurevich (1999) in the form of the Lighthill-Gueron-Liron theorem, are based on 



singular Stokeslet distributions, but employ integration to remove the near- field singularity, so that 
the local part of the integral equation is replaced by resistance coefficients. A related but different 
approach is evident in the work of Tornberg and Shelley ( 2004 ) on the modelling of complex suspensions 
of flexible fibres. These techniques have proved very successful in providing an efficient fluid mechanics 



solver as part of a model of cilia or fibre interactions with the ambient fluid ( Tornberg and Shelley 2004 



Gueron and Levit-Gurevich, 1999 2001 ). An advantage of the regularised Stokeslet method is that the 
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modelling of particles, filaments and surfaces is unified within a single mathematical framework. The 
present study is intended to expand the scope of the technique to problems which may previously have 
been computationally difficult, while retaining its ease of implementation. It should be noted that 
with sophisticated modern parallel computation, direct solution of the fluid flow and solid mechanics 



equations for large arrays of cilia is possible, see for example Mitran (2007). 



In this study we considered the flow driven by prescribed cilia and flagellar motions. In order to 
understand and interpret phenomena such as the effects of viscosity, cilia coordination and cooperation, 
axoneme ultrastructure, accessory fibres and their effect on beat pattern and frequency, it is necessary 
to couple the fluid mechanics to models of elasticity and active bending moment generation. Important 
work has already been done on this by a number of authors, including those discussed above. The 
framework presented here may assist with the associated fluid mechanics computations, particularly 
for problems with many bodies, or where there is a complex geometry to be discretised, for example 



the surface of the ciliated ventral node of the developing embryo (Hirokawa et al. , 2009; Smith et al. 



2008a), or the interaction of sperm cells with the female reproductive tract. 
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The volume flow rate due to a point-force in the vicinity of a no-slip boundary The pumping effect 
of a cilium in the x\ direction can be quantified using the volume flow rate Q, 



Q(t) 



x 3 =0 J X2 



ui(x, t) dX2dX3, 



(19) 



the choice of x\ being arbitrary due to conservation of mass. The mean volume flow rate over one 
beat cycle (1/T) Jq Q(t) dt is denoted Q. 

The volume flow rate due to a Stokeslet singularity By near to a no-slip surface was given by 



Liron (1978). Our definition of the Stokeslet differs in that we do not include the leading 1/(8717/) 
factor, so in our notation we have 



£3=0 J X2 



jBij(x,$) drc 2 dx 3 



if 
if 



1, 

2,3, 



(20) 



where x\ and £1, £2 are arbitrary. 

By conservation of mass, the volume flow rate does not depend on x\. Taking x\ and hence 
r to be large, we note that the regularisation error \S\j — Sij\ = 0(e 2 /r 3 ), and the regularisation 
error for the other images decays at least as rapidly as 0(e 2 /r 4 ). Hence the regularisation error 
\B\ - — B\j\ = 0(e 2 /r 3 ), the double integral of which tends to zero as x\ — > 00, and it cannot make any 



contribution to the volume flow rate. This proves that result (20) holds exactly for the regularised 
solution Btj, and so combining equations (19), (20) and (13) we have 



Mi M 2 ! 
^ 1=1 m =l Jo 



(s, t) ds. 



(21) 
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In i I 5.1 we report results in nondimensional units, with Q scaled with respect to cxL 3 . 

Straight line integrals of the regularised Stokeslet and their relation to slender body theory For fil- 



aments composed of straight line segments, for example the three-link swimmer (Purcell, 1977; Becker 



et al. 2003 Tarn and Hosoi 2007), equation (11) reduces to the calculation of straight line integrals. 



Similarly, provided the discretisation is sufficiently fine, curved line integrals can be approximated by 
a sum of straight line integrals. These integrals can be calculated explicitly, and utilised in a simi- 
lar manner to the singular Stokeslet-dipole integrals employed by Higdon (1979a). While this is not 
required to implement our computational method, in certain situations it is more efficient that using 
numerical quadrature, for example the results shown in p7|J7T]) and 5 5.1). As the radius of curvature 
of the relative to filament length decreases, it is necessary to subdivide the elements into multiple 
straight line segments, and the computational savings are less. For this reason, we use numerical 
integrals for the choanoflagellate example of 

For an element centred on $(s r ), having tangent £'(s r ), the approximating straight line segment 
is « (s — s r )£'(s r ) + £(s r ). The regularised Stokeslet integrals are of the form 



D(cc, s r 



-6s 



-5s 



S € {x,{s-S r )€'(s r )+£{ Sr )) ds. 



(22) 



We perform the coordinate transformation x\ = (xj — Cji s r))&ij, where is a rotation matrix. The 
first row (0ii, 0i2, ©13) is given by the tangent £'(s r ), with the second and third rows completing an 
orthogonal basis. The integral becomes, in 'local' coordinates, 



D L (x,s r 



6s 



-6s 



S e (x L ,se 1 ) ds. 



(23) 



The integral in the original coordinate system is given by Dij(x, s r ) = QkiD^ t {x L , s r ) 

Qij. The coordinate translation does not affect D since the integrand S depends on x and £ only 

through their difference. 

The definite integrals in the local coordinate system are D^(x L , s r ) = [l(j(x, s)] s J Ss , where lfj(x L , s) 
are the corresponding indefinite integrals. In concise notation (x,y,z) for (x^,X2,x^) these can be 
written, 



'22 



'33 



J 12 



J 23 



i 21 



-'32 



- n 



y 2 + z 2 + e 2 

e 2 + y 2 
y 2 + z 2 + e 2 

e 2 + z 2 
y 2 + z 2 + e 2 



! 13 



'31 



lj + 21og(s - x + r e ), 
+ log(s - x + r e ), 
+ log(s - x + r e ), 



z 

yz 



y 2 + z 2 



(24) 



with rt 



s) 2 + y 2 



+ z z + e' 



The diagonal matrix entries, occurring when x = £(s r ), require higher-order numerical quadrature 
than the other entries since the kernel is nearly singular at this point, and hence are particularly suited 



to analytic evaluation. These 'local' integrals are f_ s Ss S e (0, sei) ds, the components of which do not 
depend on x or s r . They have explicit form 



7->Local 
^11 



4 arctanh(l + (Ss/e) 



-1/2 



D 



Local 
22 
j->Local 
^12 



D 



Local 
33 
7-)Local 
^21 



2(1 + (5s/e) 2 y 1/2 + 2 arctanh(l + (ds/e) 



-1/2 



D 



Local 
13 



D 



Local 
31 



D 



Local 
23 



D 



Local 
32 



0. 



(25) 



15 



The local integrals also provide a link from the numerical implementation to classical slender body 
theory. Taking e as a proxy for slender body radius, and assuming that the length scale on which the 
force density varies, and the radius of curvature, are small in comparison with e, then it is reasonable 
to take e <C < 1. The 'local' integrals then give the dominant contribution of the Stokeslet 
distribution in the calculation of the force density, corresponding to the resistance coefficients of Gray 
and Hancock theory. Moreover, rescaling e = e/(25s), we have the asymptotic dependence on log(l/e) 



for e <C 1, as found by Batchelor (1970). 



Test problems Two test problems are considered: the translation of a fibre normal to its axis, and 
the translation of a sphere, both in an infinite fluid. The latter problem was previously considered 



extensively by Cortez et al. (2005) 



.1 Test case 1: force on an axially-translating slender body 



A classical slender body theory result (see for example Batchelor 1970, noting the slightly different 



notation) is that the magnitude of the force on a rod of length I and radius a, translating with 
velocity u parallel to its axis is, to leading order ~ 2ir fiul / log(Z / a) . Nondimensionalising with force 
scale fiul, we calculate a numerical solution using a line distribution of regularised Stokeslets, with 
the regularisation parameter chosen as e = a/ 1. We compare the two numerical implementations 
given in equations ([8]) and (10) with a/l = 0.01, and examine the force at the midpoint of the rod. 
For the latter case, integrals are calculated analytically, as described in [7j The 'boundary element' 
implementation ( 10 ) gives a relative error of less than 3.5% compared with the asymptotic estimate for 
N = 5, 10, 20, 40, 80, 160. The original implementation Q gives relative error of over 35% for N = 10, 
and it is necessary to use N = 80 degrees of freedom to give similar accuracy to the boundary element 
version. 

Benchmarking the two algorithms on a desktop PC with iV = 20 for the boundary element version 
and N = 80 for the original implementation, the former requires less than half the CPU time, using 
analytic kernel integration. If we repeat the calculation with M = 5 equally-spaced parallel rods, 
the boundary element version requires approximately one sixtieth of the CPU time compared with 
the original implementation. The reason for the significant saving is that the matrix setup cost is 
0(M 2 N 2 ), and the direct solution cost is 0(M 3 N 3 ). 



.2 Test case 2: force on a translating sphere 

To test the method with a surface distribution, we compute the drag force F = u • f s f dS^ on a 
translating sphere in an infinite fluid. Stokes' law gives this as exactly F = 6ivfj,ua, where a is the sphere 



radius. Nondimensionalising with force scale fma, we compare the boundary element method ( 10 ) with 



the results of Cortez et al. (2005). We use the same geometric discretisation, projecting the faces of 



a cube onto its surface, and subdividing these faces with a square mesh, as shown in figure [7j The 
square mesh provides a natural parametrisation of each element, which we denote (a, (3). For the 
boundary element method, numerical integration is performed using two dimensional Gauss-Legendre 
quadrature with the parametrised coordinates, taking 12 x 12 points for the 'near-singular' cases q = r 



in equation (10) and 4x4 points otherwise. The surface metric is calculated as the magnitude of 



the vector product \d£/da A d£/d/3\, as described in Pozrikidis (2002). All of the collocation and 
quadrature points lie on the curved surface of the sphere. 

For problems involving surface integrals, refining the mesh by halving the element width results 
in four times the number of degrees of freedom, 16 times the matrix setup and storage cost, and at 
least the same cost increase for an iterative linear solver. With a conventional direct solver, the cost 
increases with 0(iV 3 ), i.e. 64 times. For this reason, a method that gives accurate results with a 



relatively coarse mesh is beneficial. Consistent with Cortez et al. (2005), we compute solutions using 



the generalised minimum residual method with zero initial guess, although for the cases ./V ^ 6x12x12; 
solution by LU factorisation is possible in a few seconds on a desktop PC. 
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(a) 6 x 4 x 4 elements 
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Figure 7: Examples of computational meshes used for the test case in 7 2 table [TJ with (a) 6 x 4 x 4 
elements and (b) 6x 12 x 12 elements. The mesh shown in (a) is used in to represent the choanoflagellate 
cell body in the results shown in ^ 5.2 figures [4j [5] and [6} 



Table [jj, is based on the results with regularisation parameter e = 0.01 reported by Cortez et al. 
(2005) using equation (|8|), although we include relative percentage error and the number of scalar 



kernel evaluations necessary to construct the matrix A. In order to obtain a relative error of less than 
1% it is necessary to use a mesh with 6 x 36 x 36 patches, giving 23328 scalar degrees of freedom for 
the unknown force components, the matrix setup requiring 5.44 x 10 s kernel evaluations. 

Table [lj) shows results calculated using equation (10), with coarser meshes and fewer degrees of 
freedom. To obtain similar accuracy, a 6 x 3 x 3 mesh with 162 DOF is sufficient. While the matrix 
setup cost per DOF is greater due to the higher-order quadrature, the greatly reduced DOF results in 
the number of kernel evaluations being three orders of magnitude smaller. The matrix solution cost 
is reduced similarly. 

The effect of varying the regularisation parameter e on the predicted drag was also examined. Using 
N = 6 x 4 x 4 and constant force discretisation, the relative error for e = 0.05 was 1.4%, reducing 
to 0.63% for e = 0.01. The error remains at approximately 0.62% for e = 0.005 and 0.0025, implying 
that this is a satisfactory value — the remaining error arises from the force discretisation. Cortez et al. 



(2005) showed that the error associated with the regularisation parameter can be estimated as 0(e 2 ). 
Hence a value of e that gives accurate results for the sphere problem will in general give accurate 
results for any solid surface for which the length scales are significantly larger than e. Moreover, once 
e has been chosen sufficiently small, accurate results may be computed for any values of e that are 
smaller, provided that the kernel quadrature is computed accurately. 
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(a), Original method 



N DOF matrix entries no. kernel rel. % error 

evaluations in drag 
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6 x 36 x 36 23328 5.44 x 10 8 5.44 x 10 s 0.849 
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(b), Boundary element method 



N 


DOF 


matrix entries 


no. kernel 


rel. % error 








evaluations 


in drag 


6x3x3 
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26244 
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6x4x4 
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1.44 x 10 6 
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6x6x6 
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6x9x9 
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6 x 12 x 12 
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Table 1: Convergence of the original RSM results given in Cortez et al. (2005) and the boundary 



element implementation equation ( 10 ), for the canonical test problem of the drag force on a translating 
unit sphere. DOF is degrees of freedom. The regularisation parameter was fixed at e = 0.01 in both 
cases. The boundary element implementation used 12 x 12-point Gauss-Legendre quadrature for 
'near-singular' integrals, corresponding to q = r in equation (10), and 4 x 4-point quadrature for all 
other integrals. 
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